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1.  Introduction 


The  present  trend  of  aircraft  gas  turbine  design  has  been  characterized  by 
significant  increase  in  cycle  pressure  ratio  and  turbine  inlet  temperatures 
required  to  provide  higher  thermal  and  propulsive  efficiencies.  Also, 
increased  interest  in  engine  performance  and  fuel  econoay  has  created 
additional  emphasis  for  improving  the  efficiency  of  gas  turbine  engines.  These 
trends  accentuate  the  need  for  improvements  in  sealing  technology  and  the 
development  of  advanced  design  and  analysis  capabilities  to  reduce  gas  path 
seal  leakage,  maintain  costly  vent  leakage  to  a  minimum,  provide  better  control 
over  sophisticated  cooling  circuits,  and  prevent  high  levels  of  seal  leakage 
into  critical  aerodynamic  locations  in  the  turbine  gas  path  which  can  result  in 
a  considerable  penalty  from  thermal  and  momentum  losses.  Advanced  gas  turbine 
engine  requirements  include  a  broad  engine  power  operating  range,  which  usually 
results  in  a  wide  range  of  seal  clearance.  In  setting  the  design  clearance, 
consideration  is  given  to  transient  differential  growth,  maneuver  deflections, 
mechanical  and  thermal  growths,  eccentricity,  and  manufacturing  tolerances. 
However,  with  variable  geometry  engines  and  multiple  role  applications,  the 
engine  seals  will  not  always  operate  at  the  design  clearance  nor  provide 
minimum  leakage  across  the  operating  spectrum.  Improved  seal  design  and 
analysis  capabilities  developed  to  address  this  problem  would  have  a  major 
beneficial  impact  upon  the  design.  The  result  of  increased  cycle  pressure 
ratio  on  typical  labyrinth  seal  leakage  with  no  change  in  seal  design 
technology  is  nearly  linear,  roughly  doubling  as  the  pressure  ratio  is 
doubled.  These  increases  are  significant,  particularly  when  the  number  of 
seals  in  a  gas  turbine  are  considered.  Gas  turbines  require  a  variety  of 
labyrinth  seal  designs.  The  seal  configuration  selected  for  a  given 
application  is  based  on  the  purpose  of  the  seal  and  satisfying  design  criteria 
that  Includes  the  following  considerations:  axial  envelope  available,  axial 
travel,  clearance  range,  potential  wear,  system  sensitivity  to  seal 
clearances,  cooling  flow  requirements,  sensitivity  to  damage  in  handling, 
assembly  requirements,  and  pressure  ratio. 

Labyrinth  seals  are  used  throughout  a  gas  turbine  engine,  including: 
compressors  and  turbine  airfoil  end  seals,  bearing  compartment  seals,  and  flow 
system  seals  to  prohibit  or  control  flow.  The  purposes  of  these  seals  are  not 
always  the  same.  Labyrinth  seals  used  in  the  flow  path  are  intended  to 


minimize  end  leakage*  Bearing  coapartaent  seals  are  intended  to  keep  the  oil 
in  the  bearing  coapartaent  and  to  alniaise  the  aaount  of  leakage  and  heat 
addition  to  the  oil*  Thrust  balance  labyrinth  seals  are  located  radially  to 
provide  a  desired  off  setting  axial  load  coaponent  to  reduce  bearing  loads  to 
the  design  level.  Other  flow  system  network  seals  have  several  functions 
Including:  controlling  leakage  flows  either  to  a  aininua  or  to  a  level  to 
satisfy  disc  puaping  and  thus  prevent  hot  gaa  recirculation  in  a  cavity, 
controlling  cavity  pressure  to  reduce  axial  bearing  loads,  or  preventing 
excessive  leakage.  Seal  geometry  variables  Include  knife  edge  thickness  and 
sharpness,  clearance,  knife  pitch,  cavity  depth  and  shape,  nuaber  of  knives, 
step  height,  knife  location  on  the  step,  and  knife  angle.  Aerodynamic  u 
parameters  that  must  be  considered  in  seal  design  include  rotational  speed, 
pressure  ratio,  temperature  and  Reynolds  number. 

Seal  leakage  into  the  flow  path  has  three  loss  eleaents:  thermodynamic 
(bypassing  coabustor);  aerodynamic  (re-entry  to  flow  path  will  alter  the 
design  velocity  diagram  for  the  compressor  or  turbine  airfoils);  parasitic 
(rotational  pumping  power  required  to  bring  leakage  up  to  disc  or  mainstream 
velocity).  Improper  sealing  can  result  in  an  inadequate  supply  of  air  to 
cooled  hardware,  which  may  produce  premature  failures.  Also,  poor  sealing 
could  result  in  high  oil  consumption  and  choking  or  possibly  fires  in  the 
bearing  compartments.  In  addition,  poor  sealing  will  cause  the  engine  to 
rematch  at  lower  than  peak  component  efficiencies. 

The  benefits  of  Improved  sealing  are  wisely  recognized  and  show  the 
significant  improvements  of  reduced  seal  leakage  on  advanced  engines,  of  order 
2%  SFC  reduction  for  a  4Z  reduction  in  leakage. 

The  results  of  improving  labyrinth  seal  technology  and  reducing  leakage 
may  also  be  expressed  in  terms  of  cost  trade-offs  with  other  gas  turbine 
components.  These  trade-offs  have  shown  that  sealing  improvements  to  achieve 
the  same  level  of  compressor  or  turbine  component  efficiency  change  are  more 
economical.  In  addition,  there  is  a  significant  amount  of  improvement 
available  for  sealing  efficiency,  whereas  the  state-of-the-art  compressor  and 
turbine  component  efficiencies  have  relatively  less  room  for  improvement. 

If  a  labyrinth  seal  design  is  to  be  successful  for  the  application 
Intended,  an  accurate  seal  design  and  analysis  model  is  necessary.  The  design 
and  analysis  capabilities  available  today  rely  heavily  on  empirical 


relationships,  which  severely  limits  the  application  range.  Recently,  advanced 
numerical  techniques  based  upon  solution  of  the  governing  flow  equations,  the 
Navier-Stokes  equations,  have  been  used  for  the  labyrinth  seal  analysis  [1,2]. 
These  studies  have  been  limited  to  cases  with  two  coordinate  directions,  i.e. , 
cases  which  assume  two-dimensional  planar  symmetry  or  axial  symmetry.  For 
several  of  the  axisymmetric  cases  the  effect  of  rotation  was  considered  by  the 
imposition  of  a  rotating  boundary  and  the  additional  solution  of  the  swirl 
momentum  equation.  These  efforts  have  clearly  shown  the  ability  of  the  Navier- 
Stokes  approach  to  simulate  the  very  complex  seal  flow  field  for  a  variety  of 
practical  flow  configurations.  A  very  important  current  problem  concerns  the 
damping  performance  of  seals  in  triiich  the  gap  height  at  any  azimuthal  station 
varies  with  time.  This  can  be  a  result  of  eccentricity,  system  vibrations, 
etc.,  and  represents  a  very  challenging  problem.  The  resulting  time  varying 
shaft  loading  can  have  a  major  role  in  damping  the  vibration,  and  thereby 
changing  the  shaft  critical  speed.  The  focus  of  the  present  effort  would  be  to 
simulate  these  time-dependent  flow  fields  due  to  a  non-constant  gap  height  via 
a  solution  of  the  time -dependent  Navier-Stokes  equations. 

2.  Phase  I  -  Technical  Objective 

The  overall  objective  of  the  present  effort  was  to  develop  a  Navier- 
Stokes  analysis  which  would  be  applied  to  two-dimensional,  time-dependent 
labyrinth  configurations.  While  it  is  realized  that  the  eccentricity  or  forced 
vibration  problems  are  inherently  three-dimensional  in  nature,  it  was  felt 
that  the  potential  of  the  Navier-Stokes  approach  could  be  assessed  through 
consideration  of  a  two-dimensional  seal  geometry  in  which  the  gap  heights  vary 
in  a  prescribed  manner  with  time.  An  existing  code  with  time  invariant 
geometry  was  extended  to  apply  to  this  problem  and  calculations  were  run  to 
assess  the  procedure.  The  flow  fields  were  interrogated  to  determine  if  the 
approach  gives  the  qualitatively  expected  effect  of  gap  time-dependence.  This 
approach  was  taken  because  of  the  limited  time  and  resources  available  under 
the  Phase  I  effort.  Under  the  Phase  II  effort,  the  full  three-dimensional  time 
dependent  approach  would  be  used  to  simulate  the  eccentricity  and  forced 
vibration  problems.  Details  of  this  effort  follow. 


3.  Matheaat leal  Analysis 


The  flow  field  within  a  labyrinth  seal  is  governed  by  the  Navier-Stokes 
equations,  and  in  conjunction  with  a  suitable  turbulence  model  a  solution  of 
the  time-dependent  form  of  these  equations  would  serve  to  predict  the  flow 
field  for  both  laminar  and  turbulent  flows.  As  the  present  code  is  configured, 
the  mixing  length  mdoel  assumes  that  the  flow  is  fully  turbulent,  i.e.,  there 
is  no  transition  from  laminar  to  turbulent  flow.  The  two-equation  model, 
however,  does  have  an  inherent  transition  model  within  the  governing  equations 
and  hence,  automatically  considers  the  physics  of  transition.  The  form  of  the 
Navier-Stokes  equations  expressed  in  the  more  common  coordinate  systems  cbn  be 
found  in  standard  fluid  dynamic  texts  [3],  and  the  equations  themselves  have 
been  derived  in  general  tensor  form  by  Walkden  [4]  for  viscous  flow. 

The  form  considered  in  this  investigation  is  based  upon  density,  p,  and 
the  radial  and  axial  velocities,  u  and  w,  as  dependent  variables  and  continuity 
and  the  momentum  equations  written  in  the  cylindrical  polar  coordinates  as 
governing  equations.  In  this  approach  the  continuity  and  momentum  equations 
are  solved  in  the  form 
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and  J  is  the  Jacobian  of  the  inverse  transformation 
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Further,  the  coefficients  f$£,  Yi»  Cl  are  given  by 


8  *■  _1  ,  $2  m  _i»  ^ 

r  r 


Yx  "  l>  Y2  ■  JL.  Y3  “  1 

r 


c  -  _L  .  c2  -  I,  c3  *  1 

m  r 


and  m  ■  1  for  all  equations  except  the  x2  -  direction  momentum  equation,  for 
which  m  «  2.  The  vector  variables  used  in  Eq.  1  are  defined  as 
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Note  that  the  velocity  components  (Oj,  U2,  U3)  are  the  cylindrical-polar 
velocity  components,  and  is  the  stress  tensor  written  in  cylindricalr-polar 
coordinates.  The  molecular  and  turbulent  stress  tensors  may  be  written  £S 
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The  derivatives  required  in  Gqs.  (8)  -  (9)  oust  be  expressed  in  terms  of  the 
computational  coordinates  yj  using  the  chain  rule* 

Finally,  the  vector  £  contains  the  additional  curvature  terms  due  to 
the  cylindrical-polar  coordinate  system. 
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Since  the  present  effort  involves  the  problem  of  turbulent  flow,  a 
turbulence  model  suitable  for  this  problem  is  necessary.  Using  Favre  averag¬ 
ing  [5] ,  the  governing  equations  are  then  identical  to  the  laminar  equations 
with  velocity  and  density  being  taken  as  mean  velocity  and  viscosity  as  the 
sum  of  the  molecular  and  turbulent  viscosity.  Several  models  of  varying 
sophistication  are  available  for  the  turbulent  viscosity.  The  code  being 
used  for  the  present  effort  contains  both  a  mixing  length  turbulence  model 
and  a  two-equation  turbulence  model.  This  latter  approach  is  based  upon 
solution  of  the  turbulence  energy  and  dissipation  equations  in  conjunction 
with  the  mean  flow  equations.  In  the  present  effort,  transition  from  laminar 
to  turbulent  flow  occurs  far  upstream  of  the  first  blade  and,  in  fact,  the 
domain  where  the  computation  is  performed  is  fully  within  the  turbulent 
region.  Although  eventually  a  two-equation  model  may  be  required  for  this 
problem,  previous  experience  (Refs.  1  and  2)  has  shown  that  many  of  the 
essential  features  of  the  seal  flow  field  can  be  captured  with  a  numerical 
simulation  using  a  mixing  length  approach.  Therefore,  since  in  this  effort 
transition  was  not  of  importance  and  due  to  the  fact  that  the  use  of  a 
two-equation  model  would  increase  run  time  by  approximately  50%,  a  mixing 
length  model  was  used.  The  mixing  length,  £m,  was  computed  from  [6]. 
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where  K  is  Che  van  Karman  constant,  y  is  the  distance  normal  to  the  wall,  A  is 
the  boundary  layer  thickness  and  D  is  a  sublayer  damping  factor  which  is 
determined  from  the  analysis  of  Van  Driest  [7].  The  turbulent  viscosity,  u^, 
then  was  obtained  by 


vT  -  pv/  D:D 

where  D:D  is  the  second  invariant  of  the  mean  flow  rate  of  deformation  tensor. 
With  the  mixing  length  model  used,  the  mixing  length  will  be  zero  on  the  wall 
and  will  monotonically  increase  (not  linear)  through  the  boundary  layer 
becoming  a  maximum  of  0.09 A  at  the  edge  of  the  boundary  layer  and  remaining 
constant  thereafter.  D:D  is  a  maximum  on  the  wall  and  will  become  zero  at  the 
edge  of  a  boundary  layer.  Thus,  the  combined  effort  is  to  produce  a  turbulent 
viscosity  of  zero  on  the  wall,  increasing  to  a  maximum  with  the  boundary  layer 
and  becoming  zero  at  the  edge  of  the  boundary  layer. 

A.  Numerical  Analysis 

The  numerical  procedure  used  to  solved  the  governing  equations  is  a 
consistently  split  linearized  block  implicit  scheme  originally  developed  by 
Briley  and  McDonald  [8J  and  embodied  in  a  computer  code  termed  MINT,  an 
acronym  for  Multidimensional  Implicit  Nonlinear  Time -Dependent.  The  basic 
algorithm  has  been  further  developed  and  applied  to  both  laminar  and  turbulent 
flows.  Since  the  scheme  has  been  described  in  detail  in  several  publications 
available  in  the  open  literature,  it  will  not  be  detailed  here.  Rather,  only  a 
brief  outline  of  the  procedure  will  be  given  in  the  following. 

The  governing  equations  are  replaced  by  an  implicit  time  difference 
approximation,  optionally  a  backward  difference  or  Crank -Nicolson  scheme. 

Terms  involving  nonlinearities  at  the  implicit  time  level  are  linearized  by 
Taylor  series  expansion  about  the  solution  at  the  known  time  level,  and 
spatial  difference  approximations  are  introduced.  The  result  is  a  system  of 
multidimensional,  coupled  (but  linear)  difference  equations  for  the  dependent 
variables  at  the  unknown  or  implicit  time  level.  To  solve  these  difference 
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equations,  the  Oouglas-Gunn  procedure  for  generating  alternating-direction 
Implicit  (ADI)  splitting  schemes  as  perturbations  of  fundamental  implicit 
difference  schemes  is  introduced  in  its  natural  extension  to  systems  of  partial 
differential  equations.  This  ADI  splitting  technique  leads  to  systems  of 
coupled  linear  difference  equations  having  narrow  block-banded  matrix 
structures  which  can  be  solved  efficiently  by  standard  block-elimination 


methods. 


The  method  centers  around  the  use  of  a  formal  linearization  technique 
adapted  for  the  integration  of  initial-value  problems.  The  linearization 
technique,  which  requires  an  implicit  solution  procedure,  permits  the  solution 
of  coupled  nonlinear  equations  in  one  space  dimension  (to  the  requisite  degree 
of  accuracy)  by  a  one-step  noniterative  scheme.  Since  no  iteration  is  required 
to  compute  the  solution  for  a  single  time  step,  and  since  only  moderate  effort 
is  required  for  solution  of  the  implicit  difference  equations,  the  method  is 
computationally  efficient;  this  efficiency  is  retained  for  multidimensional 
problems  by  using  ADI  matrix  splitting  techniques.  The  method  is  also 
economical  in  terms  of  computer  storage,  in  its  present  form  requiring  only 
two  time  levels  of  storage  for  each  dependent  variable.  Furthermore,  the 
splitting  technique  reduces  multidimensional  problems  to  sequences  of 
calculations  which  are  one -dimensional  in  the  sense  that  easily-solved  narrow 
block-banded  matrices  associated  with  one-dimensional  rows  of  grid  points  are 
produced.  Consequently,  only  these  one-dimensional  problems  required  rapid 
access  storage  at  any  given  stage  of  the  solution  procedure,  and  the  remaining 
flow  variables  can  be  saved  on  auxiliary  storage  devices  if  desired.  Since 
each  one-dimensional  split  of  the  matrix  produces  a  consistent  approximation 
to  the  original  system  of  partial  differential  equations,  the  scheme  is  termed 
a  consistently  split  linearized  block  implicit  scheme.  Consistent  splitting 
has  been  shown  by  a  number  of  authors  (9]  to  considerably  simplify  the 
application  of  the  intermediate  split  boundary  conditions. 

The  present  Phase  I  effort  applied  the  LBI  scheme  to  a  time-dependent  body 
fitted  coordinate  system  simulating  flow  through  a  seal  with  a  time  varying  gap 
height.  At  every  time  step,  a  boundary  fitted  coordinate  system  was 
constructed  to  accommodate  the  time  dependent  characteristic  of  the  boundary 
(rotor).  A  simple  Eulerian  height  function  was  used  to  represent  the  rotor 
surface.  Various  techniques  to  describe  a  moving  surface  can  be  found  in  Chan 
and  Banerjee  [10].  This  function  is  defined  as  the  distance  from  a  reference 
line  and  is  prescribed  by  a  sinusoidal  function  in  time. 
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For  numerical  calculations,  the  height  function  is  discretized  into 
points,  sometimes  called  markers.  These  markers  are  then  used  to  define  the 
time  dependent  boundaries.  Since  a  boundary  fitted  coordinate  system  is  used, 
the  boundary  is  part  of  the  grid  system  at  all  times.  For  this  reason,  flow 
variables  at  the  moving  boundary  are  well  defined  and  application  of  the 
boundary  conditions  at  these  boundaries  is  straightforward. 

The  numerical  procedure  can  be  summarized  in  four  steps:  1)  the  governing 
equations  of  the  fluid  and  their  corresponding  boundary  conditions  are 
transformed  into  a  boundary  fitted  coordinate  system  according  to  the  shapes 
of  the  surface;  2)  the  governing  equations  are  solved  in  the  transformed 
coordinates  using  the  LBI  scheme;  3)  locations  of  the  surface  are  advanced 
using  the  prescribed  motion  of  the  rotor;  and  4)  advance  time  step  and  repeat 
the  sequence. 


5.  Results 

The  work  performed  in  the  present  Phase  I  program  can  be  identified  by 
three  tasks  : 

I.  Confirmation  of  the  numerical  capabilities  of  the  existing  MINT  code  to 
handle  time-dependent  sample  10-20X  scale  labyrinth  seal  geometries. 

II.  Numerical  calculations  using  realistic  flow  conditions. 

III.  Assessment  of  the  Navier-Stokes  approach  ,to  the  three-dimensional,  time 
dependent  labyrinth  seal  problem. 

These  three  tasks  are  now  discussed  in  detail. 

Task  I: 

As  the  first  task  in  the  Phase  I  effort,  the  existing  vectorized  coding 
was  verified  for  cylindrical-polar  and  time-dependent  coordinates.  These  are 
the  main  features  of  the  vectorized  code  which  have  not  yet  been  verified  but 
are  required  for  the  seal  calculations  performed  here.  The  code  has  been 
written  as  a  general  geometry  code  in  which  the  computational  coordinates  are  a 
function  of  the  Cartesian  or  polar  coordinates  and  time.  Verification 
consisted  of  two  parts.  In  the  first  part,  the  cylindrical  polar  coordinates 
option  was  validated  through  a  test  case  which  has  an  analytical  solution  for 
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the  fully  developed  laminar  pipe  flow.  The  fluid  flow  inside  a  pipe  was  used 
to  test  the  cylindrical  polar  coordinates  option  of  the  code.  For  this  case, 
the  Hagen-Poiseuille  parabolic  velocity  distribution  was  obtained.  The 
calculations  agreed  to  within  less  than  0.1%  with  the  analytical  solution  of 
the  Navier-Stokes  equations.  To  assess  that  the  time  dependent  coordinate 
option  was  working  properly,  a  series  of  detailed  printouts  including 
intermediate  results,  was  made.  These  intermediate  printouts  showed  the  code 
to  be  correctly  evaluating  the  convective -like  terms  which  are  produced  by  the 
time  movement  of  the  coordinate  system. 

Task  II:  u 

In  the  second  task,  three  demonstration  calculations  were  made.  In  these 
cases  a  geometry  of  straight  tapered  seal  with  three  knives  was  used.  The 
first  calculation  obtained  a  converged  solution  with  time-independent  seal 
boundary.  The  purpose  of  this  case  was  to  exercise  task  I  on  the  seal  geometry 
and  also  serve  as  Initial  conditions  for  cases  2  and  3.  Cases  2  and  3  involved 
calculations  of  flow  field  Inside  the  straight  tapered  seal  geometry  with  a 
prescribed  time  dependent  boundary  (Figure  1).  The  rotor  was  chosen  to  be  the 
time -dependent  surface  and  the  motion  of  it  is  sinusoidal  with  time.  In  these 
two  cases,  the  moving  surface  had  an  oscillatory  frequency  of  one.  The 
amplitude  of  oscillation  of  the  moving  surface  in  cases  2  and  3  was  0.1  and 
0.2,  respectively.  These  amplitudes  correspond  to  a  changes  of  20%  and  40%  of 
the  entire  gap  height.  . 

Since  the  overall  objective  of  the  Phase  I  effort  was  to  demonstrate  the 
ability  of  the  Navier-Stokes  solver  to  calculate  the  flow  within  the  seals  of 
an  eccentric  whirling  rotor,  a  representative  seal  geometry  with  which  the 
authors  were  well  acquainted  was  chosen.  For  all  three  cases,  fluid  flow 
calculations  have  been  performed  in  a  straight  tapered  seal  configuration  with 
three  knives  (Figure  1).  The  clearance  of  each  seal  is  chosen  as  a 
characteristic  length,  6,  of  0.1  inches  (0.245  cm).  All  other  dimensions  on 
Figure  1  are  referred  as  multiples  of  this  length.  The  stator  had  a  radius  of 
5.2  Inches  (13.2  cm).  For  the  cases  run  under  this  effort,  the  governing 
equations  consisted  of  the  transformed  strearawise  and  transverse  polar  momentum 
equations  and  the  continuity  equation.  The  initial  conditions  and  boundary 
conditions  can  be  summarized  as  the  following.  The  flow  is  initially  assumed 
to  be  stagnant  and  the  back  pressure  lowered  to  the  desired  level.  The 


boundary  conditions  use  the  no-slip  conditions  on  the  walls,  the  two-layer 
model  at  the  inlet,  and  specified  static  pressure  on  the  exit  boundary.  The 
two  layer  model  divides  the  upstream  inlet  boundary  into  two  regions  (1)  a 
central  core  in  which  the  total  pressure  is  specified  and  (2)  a  boundary  layer 
region  in  which  the  pressure  is  assumed  constant  and  the  form,  but  not  the 
magnitude,  of  the  boundary  layer  velocity  profile  is  assumed.  The  assumed 
profile  is  normally  determined  by  specifying  a  distance  from  a  leading  edge  and 
a  boundary  layer  growth  rate.  Thus  a  boundary  layer  thickness  is  determined  at 
the  computational  inlet  plane  and  a  boundary  layer  velocity  profile  is  either 
user  specified  or  determined  by  a  variety  of  analytical  or  semi -analytical 
techniques  (e.g.,  (11-12]).  In  this  study,  the  method  of  Maise  and  McDoAald 
[11]  was  used.  In  essence,  the  boundary  layer  thickness  and  the  form  of  the 
velocity  profile  represent  the  upstream  history  of  the  flow  before  entering  the 
region  of  computation.  The  edge  velocity  was  determined  from  the  specified 
inlet  total  pressure  and  the  static  pressure  determined  via  the  time-evolving 
solution.  For  the  cases  considered,  the  pressure  ratio  was  2.0.  The  inlet 
boundary  layer  profiles  on  the  rotor  and  land  entrance  were  assumed  to  have  a 
thickness  of  0.05  inches  (0.0122  cm)  and  a  skin  friction  coefficient  of  0.005. 
The  turbulence  viscosity  was  assumed  to  be  modeled  by  a  mixing  length  model 
previously  discussed.  Once  a  steady  state  was  achieved,  the  rotor  is  allowed 
to  move  in  the  radial  direction. 

Figure  2  and  4  show  the  velocity  vector  plots  and  pressure  contours  inside 
the  blade  gaps  for  case  1.  Figure  3  shows  the  axial  velocity  profile  along  the 
first  coordinate  line  in  the  first  blade  gap.  In  this  case,  the  flow  remained 
subsonic  except  in  a  small  axial  region  (656  to  676)  near  the  third  knife 
blade,  where  6  refers  to  the  gap  shown  in  Fig.  1.  The  maximum  Mach  number  was 
1.03  at  an  axial  location,  656,  in  the  third  blade  gap.  This  was  approximately 
where  the  dimensionless  pressure  is  0.48.  A  large  recirculation  exists 
downstream  of  the  last  blade  and  recirculating  zones  exist  in  the  two  cavity 
regions,  as  was  expected.  A  small  recirculating  vortex  was  found  at  the  first 
blade  surface.  Along  all  blade  surfaces,  low  pressure  points  were  found  at  the 
leading  edge  of  each  blade.  In  the  third  blade  surface,  a  low  pressure  point 
was  also  observed  at  the  trailing  edge.  At  the  exit  area  of  the  first,  second 
and  third  blade,  the  pressures  are  0.64,  0.60  and  0.48  of  the  inlet  pressure, 
respectively. 
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In  cases  2  and  3,  the  radial  aotion  of  the  rotor  was  assured  to  have  an 
oscillatory  frequency  (in  tine)  of  unity.  To  study  the  effects  due  to  this 
moving  boundary,  pressure  and  axial  velocity  near  the  middle  of  the  third 
knife  blade  surface  were  traced  in  time  both  on  the  blade  and  the  rotor  at  the 
position  shown  in  Fig.  1.  The  results  are  plotted  in  Figures  5  and  19  for 
these  cases.  Figure  5  shows  the  pressure  and  axial  velocity  for  case  2.  In 
this  case,  the  amplitude  of  oscillation  for  the  rotor  was  0.1.  This 
corresponds  to  a  gap  height  change  of  0.96  to  1.16.  Approximate  periodic 
solutions  were  obtained  after  five  transient  cycles.  In  this  case  and  the  next 
set  of  calculations,  the  radial  velocity  of  the  fluid  on  the  rotor  was  imposed 
by  the  prescribed  motion  of  the  rotor.  The  amplitude  of  the  oscillatory  *-rotor 
was  0.1.  For  this  location,  the  pressure  variation  at  the  rotor  surface  was 
0.40  of  that  at  the  third  blade  surface  and  these  two  pressures  had  a  phase 
shift  of  about  45°.  Figures  6-9  show  the  coordinates  in  the  three  blade  gaps 
at  quarter  cycles  of  the  motion  of  the  rotor.  Notice  that  the  moving 
coordinates  occur  in  the  gap  regions  only;  this  is  to  reduce  any  unnecessary 
numerical  disturbance  in  other  regions.  Figures  10  -  13  show  the  velocity 
vectors  in  the  blade  gaps  at  different  quarter  cycles  of  the  movement  of  the 
rotor.  Figure  14  shows  the  axial  velocity  profiles  along  the  first  coordinate 
line  in  the  first  blade  gap  at  different  quarter  cycles  of  the  movement  of  the 
rotor.  Small  recirculating  vortices  were  found  at  the  first  stator  blade 
surface  in  all  situations.  The  entire  flow  remained  subsonic  except  in  a  small 
region  near  the  exit  of  the  third  stator  blade.  The  Mach  number  in  this  region 
varied  from  1.06  to  1.1,  with  the  largest  value  occurring  when  the  rotor  was  at 
its  first  quarter  cycle  (i.e.,  when  the  gap  height  is  one  and  is  on  the  way  to 
a  smaller  height).  This  maximum  Mach  number  occurred  at  an  axial  distance  of 
666.  Figures  15  -  18  show  the  pressure  contours  in  the  blade  gaps  at  different 
quarter  cycles  of  the  movement  of  the  rotor.  The  time  dependent  nature  of  the 
flow  field  is  clearly  evident  in  these  figures.  In  particular.  Figures  15  and 
17  show  the  same  instantaneous  gap  height  with  the  gap  increasing  in  Figure  17 
and  decreasing  in  Figure  15.  The  pressure  contours  for  these  times  show 
significant  variation. 

In  the  last  case  (case  3),  the  amplitude  of  oscillatory  rotor  was  chosen 
to  be  0.2.  This  corresponds  to  a  change  of  gap  height  from  0.86  to  1.26.  The 
pressure  and  axial  velocity  variation  at  the  starred  location  is  shown  in 
Figure  19.  Figures  20  -  23  show  the  velocity  vectors  in  the  blade  gaps  at 
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different  quarter  cycles  of  the  movement  of  the  rotor.  Figure  24  shows  the 
sxlal  velocity  profiles  along  the  first  coordinate  line  in  the  first  blade  gap 
at  different  quarter  cycles  of  the  movement  of  the  rotor.  Small  recirculating 
vortexes  were  found  at  the  first  stator  blade  surface  In  the  third  and  fourth 
quarter  cycles  only.  The  Mach  nuaber  varies  froa  1.10  to  1.13,  with  the 
largest  value  at  the  third  quarter  cycle  (when  the  gap  height  is  one  and  is  on 
the  way  to  larger  gap  height).  This  aaxiaua  Mach  nuaber  occurred  at  an  axial 
distance  of  676.  Figures  25  -  28  show  the  pressure  contours  in  the  blade  gaps 
at  different  quarter  cycles  of  the  aoveaent  of  the  rotor. 

In  summary,  of  the  calculations  discussed,  the  first  was  for  a  steady 
state  solution.  This  was  also  used  as  an  initial  condition  for  the  subsequent 
time-dependent  runs.  The  numerical  technique  utilised  in  this  first  case* was 
to  obtain  a  steady  state  solution  in  the  most  computationally  efficient  manner 
without  regard  for  transient  accuracy.  This  was  accomplished  by  a  matrix 
preconditioning  technique  [13]  which  can  be  viewed  as  taking  a  fictitious 
variable  time  step  through  the  spatial  field.  For  cases  2  and  3,  which 
involved  time-dependent  solutions,  the  physical  time  steps  were  used  to 
maintain  transient  accuracy.  A  constant  time  step  of  0.005  was  chosen  for 
these  cases.  The  periodic  solution  obtained  indicates  this  time  step  choice  to 
be  a  viable  one. 

The  converged  solution  for  the  first  case  was  obtained  in  300  iterations. 
The  pressure  was  lowered  to  the  ambient  value  equal  to  50  percent  of  the 
stagnation  value  over  25  time  steps,  and  then  the  case  was  run  to  the  steady 
state  solution.  For  cases  2  and  3,  with  a  oscillating  land  of  unit  frequency, 
200  time  steps  were  required  to  complete  one  oscillatory  cycle.  The  CPU  time 
for  each  time  step  was  2.25  seconds  with  75*150  grid  points.  A  fully 
vectorized  version  of  MINT  code  was  used  for  the  present  calculations. 

Task  III; 

The  objective  of  the  present  Phase  I  effort  was  to  perform  a  preliminary 
assessment  of  the  capability  of  a  Navier-Stokes  code  to  simulate  flow  physics 
of  a  seal  with  varying  gap  and  to  assess  the  practicality  of  such  a  procedure. 
Considering  the  latter  item  first,  a  periodic  solution  was  obtained  in  1125 
seconds  of  CRAY-XMP  CPU  time.  This  run  time  is  within  range  to  allow 
simulations  of  large  numbers  of  selected  cases.  Further,  since  no  study  of 
maximum  allowable  time  step  has  been  made,  it  is  possible  that  the  allowable 
time  step  could  be  increased  considerably,  thus  shortening  required  run  time. 


* 


V 


g 

1 


-14- 


8 


wnnuiuiMJiL 


In  regard  Co  the  results,  periodic  solutions  were  clearly  obtained.  This 
is  most  dramatically  shown  by  comparison  of  pressure  distributions  for  the  same 
gap  height  at  times  180°  apart  in  the  cycle  (Figures  15  and  17).  In  regard  to 
an  approximate  assessment  of  the  results,  it  is  possible  to  consider  the 
present  two-dimensional  calculation  as  simulating  some  features  of  the  actual 
three-dimensional  case.  This  can  be  done  by  considering  the  instantaneous 
solution  at  any  time  to  correspond  to  the  solution  for  the  three-dimensional 
problem  at  the  azimuthal  location  having  the  same  gap  height.  The  net  force 
acting  on  the  rotor  can  then  be  calculated  from  the  pressure  at  the  rotor  in 
one  cycle  of  oscillation. 

The  net  forces  acting  on  the  rotor  are  calculated  for  each  case  (2  and  3). 
In  case  2,  the  calculated  net  force  (average  dynamic  centering  force)  is  3.2  N 
and  is  acting  at  an  angle  of  141.5°  (see  Figure  29).  By  assuming  an  average 
radius  of  whirl  orbit  of  0.16  inches,  the  radial  stiffness,  Kgjj-F/r,  is  12.6 
KN/m.  The  radial  stiffness  in  a  similar  situation  [14]  is  obtained 
experimentally  in  the  range  of  0  to  25  KN/m  over  a  range  of  frequencies  and 
pressure  gradients.  Our  two-dimensional  approximation  to  a  whirling  rotor  case 
thus  provides  a  qualitatively  satisfactory  result. 

6.  Conclusions 

Under  the  present  effort,  a  transient  capability  has  been  developed  which 
enables  a  designer  to  calculate  the  flow  in  a  labyrinth  seal  configuration 
with  an  eccentric  whirling  rotor.  This  has  been  done  by  solving  the  transient 
Navier-Stok.es  equations  in  a  coordinate  system  which  reflects  the  seal  gap 
variation  as  a  function  of  time.  This  new  capability  has  been  incorporated 
within  a  very  general  two-  and  three-dimensional  Navier-Stokes  procedure  which 
has  been  applied  previously  to  a  wide  variety  of  turbomachinery  flow  problems. 
Typical  examples  have  been  seals  with  and  without  swirl,  disk  pumping  cavities, 
cascades,  etc.  Although  the  analyses  have  been  confined  to  equations  written 
in  an  inertial  frame,  extension  to  a  rotating  frame  is  straighf orward . 
Therefore,  the  basic  code  can  now  effectively  simulate  a  large  number  of 
relevent  flow  fields. 

In  regard  to  the  work  considered  specifically  under  this  effort, 
calculations  were  made  in  an  axisymmetric  simulation  for  a  sinusoidally  varying 
gap  height  at  two  frequencies.  Periodic  flow  was  obtained  in  both  cases  with 
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five  cycles  of  aotion.  The  computed  flow  fields  were  examined  and  the  results 
appeared  physically  realistic.  Obviously,  a  detailed  comparison  with 
experiment  for  three-dimensional  configurations  of  interest  is  required  before 
a  quantitative  assessment  can  be  made.  The  basic  procedure  has  been  applied  to 
a  variety  of  three-dimensional  problems  and  application  to  a  seal  problem  which 
has  time-dependent  gap  height  and  is  not  axlsymmetric,  l.e.,  fully  three- 
dimensional,  should  be  relatively  straightforward. 

The  present  SRA  Navier-Stokes  code  represents  a  state-of-the-art  procedure 
in  simulating  highly  resolved  complex  flow  fields.  The  procedure  is  based  upon 
an  efficient  ADI  approach  which  allows  high  near  wall  resolution  and 
significant  grid  stretching;  both  of  these  properties  are  required  to'  compute 
the  complex  flow  present  in  turbomachinery.  The  code  itself  is  highly 
vectorised,  leading  to  very  short  run  times  per  grid  point  per  time  step.  For 
steady  state  solutions,  the  code  contains  matrix  reconditioning  techniques 
which  lead  to  very  rapid  convergence.  Typically,  convergence  is  obtained 
within  100  to  300  time  steps  for  flows  having  near  wall  resolution  and  a  large 
number  of  grid  points.  These  properties  make  the  code  very  suitable  for 
application  to  complex  turbomachinery  flow  problems. 
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Pressure  and  axial  velocity  near  the  middle 
of  the  third  blade  vs.  time. 
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Figure  11.  Case  2  -  Velocity  vectors  inside  the  blade  gaps  at  t  -  9.25. 
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Figure  22.  Case  3  -  Velocity  vectors  inside  the  blade  gaps  at  t  -  4.75. 
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